Understanding the model used by lme4

lme4’s model

Simplified from Bates et al. (2015).

Linear model conditional on random effects $b$:
$$(y|b)\sim \mathcal{N}(X\beta +Zb,{\sigma }^2{I}_n)$$
$X$ is an $n\times p$ design matrix.

Distribution of $b$:
$$b\sim \mathcal{N}(0,\mathrm{\Sigma })$$
$Z$ is an $n\times q$ design matrix, $\mathrm{\Sigma }$ is positive definite.

More precisely, there can be several random effects terms $i$, $Z_i{b}_i$, and for each $Z_i$ contains a $q_i$ regressors within $r_i$ groups. The coefficients for the regressors within each group have a $q_i\times q_i$
covariance ${\mathrm{\Sigma }}_i$ across groups. I.e.
$$b_i\sim \mathcal{N}(0,I_{r_i}\otimes {\mathrm{\Sigma }}_i)$$

The example data

Load sleepstudy and restrict to 3 days and 3 subjects:

library(lme4, quietly=TRUE)
data <- sleepstudy[(sleepstudy$Days <= 2) & (sleepstudy$Subject %in% c(308, 309, 310)), ]
data

	Reaction	Days	Subject
1	249.5600	0	308
2	258.7047	1	308
3	250.8006	2	308
11	222.7339	0	309
12	205.2658	1	309
13	202.9778	2	309
21	199.0539	0	310
22	194.3322	1	310
23	234.3200	2	310

Model, variables are (Intercept) and Days, both fixed and random:

model <- lmer(Reaction ~ Days + (Days | Subject), data)
summary(model)

Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (Days | Subject)
   Data: data

REML criterion at convergence: 65.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3294 -0.3361 -0.0463  0.5375  0.9859 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 Subject  (Intercept) 808.5    28.43         
          Days        122.8    11.08    -0.56
 Residual             140.0    11.83         
Number of obs: 9, groups:  Subject, 3

Fixed effects:
            Estimate Std. Error t value
(Intercept)  221.403     17.561  12.607
Days           2.792      8.016   0.348

Correlation of Fixed Effects:
     (Intr)
Days -0.585

We now try to extract the parts of the model in the equation above (constants and estimates).

fixed effects

model matrix

getME(model, "X")

   (Intercept) Days
1            1    0
2            1    1
3            1    2
11           1    0
12           1    1
13           1    2
21           1    0
22           1    1
23           1    2
attr(,"assign")
[1] 0 1
attr(,"msgScaleX")
character(0)

c(getME(model, "n"), getME(model, "p"))

[1] 9 2

parameter estimates

# getME(model, "beta")
getME(model, "fixef")

(Intercept)        Days 
 221.402556    2.791767 

random effects

model matrix

getME(model, "Z")

9 x 6 sparse Matrix of class "dgCMatrix"
   308 308 309 309 310 310
1    1   .   .   .   .   .
2    1   1   .   .   .   .
3    1   2   .   .   .   .
11   .   .   1   .   .   .
12   .   .   1   1   .   .
13   .   .   1   2   .   .
21   .   .   .   .   1   .
22   .   .   .   .   1   1
23   .   .   .   .   1   2

c(getME(model, "n"), getME(model, "q"))

[1] 9 6

residual standard deviation

# getME(model, "sigma")
sigma(model)

[1] 11.83223

getME(model, "theta")

     Subject.(Intercept) Subject.Days.(Intercept)             Subject.Days 
               2.4031809               -0.5266796                0.7742044